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A limit on eccentricity growth from global 3-D simulations of 
disc-planet interactions 
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ABSTRACT 

We present high resolution 3-D simulations of the planet-disc interaction using smoothed par- 
ticle hydrodynamics, to investigate the possibility of driving eccentricity growth by this mech- 
anism. For models with a given disc viscosity (a = 0.01), we find that for small planet masses 
(a few Jupiter masses) and canonical surface densities, no significant eccentricity growth is 
seen over the duration of our simulations. This contrasts with the limiting case of large planet 
mass (over twenty Jupiter masses) and extremely high surface densities, where we find eccen- 
tricity growth in agreement with previously published results. We identify the cause of this 
as being a threshold surface density for a given planet mass below which eccentricity growth 
cannot be excited by this method. Further, the radial profile of the disc surface density is 
found to have a stronger effect on eccentricity growth than previously acknowledged, imply- 
ing that care must be taken when contrasting results from different disc models. We discuss 
the implication of this result for real planets embedded in gaseous discs, and suggest that the 
disc-planet interaction does not contribute significantly to observed exoplanet eccentricities. 
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1 INTRODUCTION 

Modern theories of planet formation all posit that planets are the 
end result of the evolution of circumstellar discs of gas and dust, 
termed protoplanetary discs. This basic idea is consistent with both 
the observation of such di s cs around young, solar-like s tars (e.g., 
ISargent & Beckwithl fl987l : lHaisch. Lada. & Lad3l200ll) . and the 
near co-planar, near circular orbits of the planets in the Solar Sys- 
tem. However, observations of extra-solar planets made over the 
last 15 years have revealed that many planetary systems do not 
share the neat structure of our own. In particular, the abundance of 
planets with eccentric orbits was unexpected given this formation 
mechanism. The observed distribution of eccentricities is nearly 
uniform between eccent ricity e = and ~ 0.6, and stretch es all 
the way up to e - 1 (e.g. lWright et alJl201ll ; lKane et alj|2012h . 

At small semimaj or-axis (a < 0. 1 au) there is a clear preference 
towards low eccentricit y (e < 0.1) orbits, and this is well explained 
by tidal circularisation dRasio etal.lll996t) . Similarly, the dearth of 
planets in high eccentricity orbits (e > 0.5) with small semi-major 
axes (< lau) is readily understood, as these planets have periastron 
distances ~ 1 R Q and consequently pass too close to their stars to 
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be long-lived. At larger semi-major axes, however, the picture is 
much less clear, and a number of scenarios have been proposed to 
to explain the apparent discrepancy between theory and observa- 
tion. These fall into three broad categories: dynamical interactions 
of planets in multi-planet systems; secular interactions with com- 
panion stars; and tidal interactions between planets and their parent 
protoplanetary discs. 

Large regions of the observed eccentricity distribution can 
be populated by invoking interactions between multiple planets, 
be it via direct close e ncount er s dFord. Havlickova. & Rasid 
1200 ll : ljuric & Tremaind 120081 ; IChatteriee et alj l2008h or 
through mean-motion resonances over longer time-scales 
dChiang. Fischer. & Thommesl l2002h . While simulations of such 
encounters are able to reproduce the observed distribution to 
a reasonable accuracy, it is unclear if such close interactions 
are frequent enough in nature to provide a universal source of 
planetary eccentricity. Another method for growing eccentricity is 
through secular interactions with inclined companion stars, which 
lead to a resonant exchange of angular momentum betwe en the 
planets and the external body jKozail 1 19621 ; llidovi 1 19621) . This 
results in long-period changes in inclination and eccentricity 
and although this mechanism seems inviting as an alternative 
explanation for the observed eccentricity distribution, numerical 
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work has s hown that it does not p roduce the correct eccentricity 
distribtion dTakeda & Rasiol l2005h . although recent work has 
found that it may explain s ome eccentric misaligned Hot Jupiters 
dNaoz. Farr. & Rasioll 20121) . Similar interactions between planets 
in the same system have also been sugge sted as a chaotic forma tion 
mechanism for highly eccentric planets dWu & Lithwickl20llh . 

The gravitational interaction between a young embedded 
planet and its parent gas disc has also been suggested as a mecha- 
nism for driving eccentricity growth. For companions with masses 
comparable to the central body it has long been known that tidal 
interactions with the disc lead to eccentri city excitation. This resul t 
has applications to stellar binaries (e.g.. lArtvmowicz et al.lll99lh 
and binary super-massive black holes (e.g. lCuadra et alj|2009l) . but 
how it extends to the more extreme mass ratios of star-planet sys- 
tems is still not clear. 

Analytical treatment of this problem begins with the consid- 
eration of resonant torques between a planet and its parent disc. 
These torques come in two main flavours, Lindblad and corotation. 
Lindblad torques occur where the epicyclic frequency in the disc is 
equal to plus or minus the frequency of a component of the planet's 
potential. Corotation resonances occur at radii where the compo- 
nent's orbital frequency matches that of the disc. For a planet with 
a circular orbit, the potential can be described by a series of com- 
ponents with the same angular velocity as the planet. This results 
in a coorbital corotation resonance with combs of Lindblad reso- 
nances to either side. When the planet's orbit is eccentric, this is no 
longer the case. A more complex picture emerges where not all res- 
onances have the sam e angular velocity, resulting in non-coorbital 
corotation resonances dGoldreich & Tremainell98rj|) . The interplay 
between these resonances, the planet that excites them, and how 
they are affected by various disc p aram eters have been di s cusse d in 
detail bv lQgilvie & Lubowl d2003l) and lGoldreich & Saril d2003l) . A 
great deal of nuance between competing effects is revealed, and in 
particular the possibility that corotation torques may begin to sat- 
urate and weaken once an initial eccentricity is attained is useful 
when considering the results o f numerical work dOgilvie & Lubowl 
l2003tlMasset & Ogilvidl2004h . However, the complex, non-linear 
nature of the planet-disc interaction makes the general problem an- 
alytically intractable, and numerical solution is required to gain a 
full understanding of eccentricity growth. 

Semi-analytic calculations combining prescriptions from 
these analytic wo r ks ha ve been somewhat inconclusive. 
iMoorhead & Adams] d2008h found eccentric damping rather 
than growth in most cases, although in the cases where they did 
find growth it was extremely strong, leading to e ~ 1 after only 
a few thousand orbits. However, as such highly eccentric planets 
would be unable to maintain an equally eccentric gap their orbits 
wou ld be circularised as they interact with coorbital disc material 
(e.g. lBitsch&KlevkOlOh . 

By contrast, numerical simulations looking specifically at 
eccentricity gro wth have so far shown mor e cons istent and 
positive results. IPapaloizou. Nelson. & Massetl J200 lL hereafter 
IPNMOll) found that relatively massive embedded companions (in 
the brown dwarf regime) undergo eccentricity growth. Lower- 
mass planets were not found to experience this growth, although 
it seems lik ely that this was for n umerical rather than physi- 
cal reasons dMasset & Ogilvid |2004) . Later simulations have in- 
deed fo und eccentricity growth, albeit at modest levels, down 
to MT,,p dD'Angelo. Lubow. & Batd 1 20061) . Extensiv e analys is of 
the behaviour and m orphology of the disc by IPNMOll and 
iKlev & Dirkser] J2006h attributes this eccentricity excitation to an 
instability launched at the 3:1 outer Lindblad resonance, which 



drives a large eccentricity at the in ner edge of the disc. For the large 
companion masses considered bv lPNMOll a wide gap is opened in 
the disc, so coorbital co-rotation resonances are not present, and 
non-coorb ital ones only operate o nce the planet's orbit is already 
eccentric. IKlev & Dirkser] J2006h also extend this analysis down 
to planets of a few M lup , and find that this mechanism still oper- 
ates down to planets of mass ~ 3 Mj up , although the magnitude of 
the eccentricity induced depends strongly on the disc viscosity and 
temperature. 

To date the majority of the numerical simulations of this prob- 
lem have been performed in only two dimensions (2-D), and all 
have used Eulerian (grid-based) methods. However, it has been sug- 
gested that a full three-di mensional (3-D) treatment weake ns the 
effect of resonant torques dTanaka. Takeuchi. & Wardll2002f) . Each 
study has also typically only considered a single disc model, with 
little consistency in the choice of parameters, and it can be eas- 
ily seen that the properties of the disc itself plays a large role in 
the evolution of a planet embedded within it. Moreover, Eulerian 
methods are not always ideal for following the dynamics of gas on 
non-circular orbits; in general, one expects Lagrangian methods to 
track eccentric orbits with greater accuracy. 

In this paper we present results of high-resolution 3-D 
smoothed particle hydrodynamics (SPH) simulations of eccentric- 
ity growth due to planet-disc interactions. In section|2]we describe 
our numerical method and the initial conditions used. In section [3] 
w e present results of our simulations, comparing them to the results 
of IPNMfjll especially, and describing the effect of various model 
parameters. Section [4] contains a discussion of both the numerical 
limitations of our simulations, and the physical interpretation of the 
results. Our summary and concluding remarks are in section[5] 



2 SIMULATIONS 

We have performed a suite of three-dimensional simulations of a 
planet embedded in a protoplanetary accretion disc. W e used a 
modi fied version of the publicly available Gadget-2 code dSpringell 
l2005h . a hybrid SPH/N-body code. SPH methods are highly suited 
to tracing the evolution of dynamical systems, as the Lagrangian 
formulation means that important orbital parameters ( energy, linea r 
& angular momentum) are naturally conserved (e.g., |Pricdl2012h . 
We treat the star and planet as point masses, and model the gas 
disc using SPH. The basic parameters of the code are unchanged 
throughout our simulations, and were set as follows. The number 
of nearest neighbours for kernel integration is 50, and smoothing 
lengths are adjusted accordingly when the number of neighbours 
changes by more than ±2. The Courant parameter used to deter- 
mine the time-step for SPH particles was set to 0.1. The gravita- 
tional softening length for the point mass particles was in each case 
set to be the same as the sink radius for the planet particle (de- 
scribed below). SPH particles are therefore accreted by the sink 
particles before they encounter gravitational softening, so this soft- 
ening has no effect on our simulations. 

We have altered the public Gadget-2 code in a number of ways 
to make it more suitable for simulating planets embedded in discs. 
Firstly, we note that the standard Barnes-Hut tree is insufficiently 
accurate when computing the gravitational forces on the N-body 
particles ( particularly the 'pla net'). Consequently, following the 
method of lcuadra et all d2009h . we removed the N-body particles 
from the tree, and instead computed their gravitational forces by di- 
rect summation. As this is done for only two particles the additional 
computational cost is not large, and ensures high accuracy for the 
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orbit of the planet. Additionally, following lCuadra et al, I l2006h . wc 
treat each N-body particle as a 'sink' for SPH particles, accreting 
the mass and momentum of the swallowed particles on to the target 
sink, ensuring conservation of linear and angular momentum to a 
high accuracy. For the planet particle, the radius within which the 
sink will accrete SPH particles is set to be 0.4 of the Hill radius, 
defined as 



^Hill - a 



Mp 
3M, 



1/3 



(1) 



for a planet and star of mass M p and M* respectively with semi- 
major axis a. The sink radius for the star particle was set to be 7/8 
of the initial inner disc radius. We have also substituted the energy 
and entropy equations in Gadget-2 with a locally isothermal equa- 
tion of state, where the sound speed c s in the disc is a function only 
of the cylindrical radius R. We choose a power-law profile for the 
disc temperatur e T(R) oc R~ i/2 , cons istent with both a linear vis- 
cosity law ( e.g., Hartmann et all 19981) a nd a standard 'flaring disc' 
model (e.g- lKenvon & Hartmann! 19871) . 



2.1 Artificial viscosity and angular momentum transport 

In order to smooth out the discontinuities associated with shocks in 
SPH, it is nec essary to make use of an arti ficial viscosity. We use 
the method of iMorris & Monaghanl i ll 9971) to dissipate energy in 
shocks and prevent particle penetration. Thi s contributes a further 
acceleration term to the equations of motion (Springel 20050 



dit 



au„ v~i — 
-2- = -^_ i m b n ab V a W a 



(2) 



W a b is the SPH kernel used by Gadget-2, a cubic spline that goes to 
zero at one smoothing length h, while Wab is the arithmetic mean 
of W ab (h a ) and Wba(h b ^ Yl a/ , describes the strength of the artificial 
viscosity and is given by 



n =1 ( -~ a " b L ' Sal > t * ab + P^cJIPab for u «b ■ r ab < 0, 
ab \ otherwise. 

with p given by 

hob u ab ' r ab 



Pab 



\r ab \ 2 + sh„ 



(3) 



(4) 



Here c7 oi , p ab and h ab represent the arithmetic means of the SPH 
sound speed, density and smoothing lengths respectively, between 
particle a and b while u a b and r a b are the differences between the 
velocity and position vectors of those particles, e = 0.0001 is a 
numerical 'safety factor', which prevents the artificial viscosity di- 
verging for very small particle separatio ns. 

In the lMorris & Mona ghan ( 199% scheme, /3 in equation [3] is 
given by p = 2a ab , an d a ab is the mean value (a a + a b )/2. We have 
used the IPricd (2004) version of this method, where a a is evolved 
according to 



(5) 



1 In this system of equations, i, j and k are vectors using summation nota- 
tion, and a and b are particle indices. 

2 Note that this differs from the 'standard' SPH definition of the smoothing 
length, where the kernel goes to zero at 2h. This difference in notation has 
no practical effect, but must be considered when comparing these equations 
to those in the SPH literature. 



(6) 



(7) 



and a m i n and a max are imposed minima and maxima, set as a global 
parameters. The source term S is given by 

S a = max(-V • «„, 0} 

and the decay time-scale is 

K 

T a = . 

where / is a dimensionless scaling parameter. We follow IPricd 
d2004l) and chose values / = 0.1, a mi „ = 0.01 and a ma x = 2.0. 
We also make use of the 'Balsara switch' (Balsara 19951) to limit 
the spurious detection of shear flow as shocks. 

This artificial viscosity prescription is necessary to handle 
shocks correctly, but we also wish to include an explicit physical 
viscosity to drive angular momentum transport in the disc. To this 
ef fect we have m odified the SPH equation of motion (equation 7 
in ISpringelll2005l) to include a Navie r-Stokes shear viscosi ty term 
(with the bulk coefficient set to zero; lLodato & Pricel2010l) : 



<K 

dt 



2 



nil, 



(8) 



Here A n is dimensionless and accounts for smoothing-length gra- 
dients: 



A, 



1 + 



3p a dh a 
and S'J is the stress tensor 



(2 \ 9«* 



8" + 



(H Qui) 



(9) 



(10) 



P and p are the SPH-estimated pressure and density respec- 
tively. The standard va riable-smoothing-length operator is given by 
dLodato & Pricel2010l) 



H 1 v / , j \ dW ab (h a ) 



(11) 



j] a is the shear viscosity, which is parametrized in terms of kine- 
matic shear viscosity v a by 



rial Pa- 



iU) 



In practice we set the viscosity such that each particle a has a 
viscosity v a which depends on only its (cylindrical) radius, which 
gives a straightforward means of implementing the power-law vis- 
cosity prescriptions described in section [2~2l 



2.1.1 Viscous ring-spreading tests 

To test this set-up, we conducted test simulations which model 
the viscous spreading of a gas ring around a point mass. In these 
tests a thin ring of initially uniform surface density is allowed to 
evolve under the acti on of viscous torques. We expect the ring 
to spread (see, e.g., IPringlel [l98lh , and by comparing our re- 
sults against those from a one-d imensional explicit scheme (e.g., 
IPringle. Verbunt. & WadeHl98fj) we can test the accuracy of our 
SPH viscosity prescription. We modelled the ring using 10 5 and 
10 6 SPH particles and four different levels of viscosity (see table 
[TJ. The ring was set up orbiting a single point mass, and was al- 
lowed to evolve for 200 orbits. 

The viscous ti me-scale for a s preading ring with constant vis- 
cosity is is given by jPringlell98ll) 
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Figure 1. Radial surface density evolution of a viscously spreading ring. Panel (a) shows a ring with 10 particles in its initial configuration (solid black line), 
after 4.9 orbits (dashed black line) and after 39.6 orbits (dash-dotted black line). The corresponding red lines show the best-fitting profile from an explicit 1-D 
ring code, plotted as a fraction of the viscous spreading time. Panel (b) is as (a), but for a ring with 10 6 particles, in its initial configuration (solid line), after 
12 orbits (dashed line) and after 90 orbits (dash-dotted line). Our Navier-Stokes viscosity approximates a uniform viscosity well except at very low resolution, 
where the 'best fit' to the 1-D ring is rather poor in both cases. The initial configuration in panel (b) is representative of the resolution we obtain in our full 
disc models described in section|2.2| 



We measure the effective viscosity in the SPH simulations by fit- 
ting the resulting surface density profiles T.(R, t) to those obtained 
from the 1-D explicit code. We perform a least-squares fit for Y.(R) 
at several different times in each simulation, and measure the effec- 
tive viscous time-scale t y . The ring spreading is initially linear with 
the 1-D model as expected, but later the approximation of constant 
viscosity becomes invalid as the artificial viscosity becomes strong 
at the ring edges. We only fit spreading times during this initial 
phase where the linear relationship exists. We are thus able to com- 
pare the imposed viscosity v- m with the measured rate of viscous 
angular momentum transport v ou t> in order to determine the accu- 
racy of our numerical s cheme. We further parametr ize the viscosity 
in terms of an effective Ishakura & Sunvaevl ( fl973h a parameter by 
assuming that v = ac^Ho at Rq. Our measured values are given in 
tableQ] and typical best fits are shown in figure[TJ(for the runs with 
v- m = 10 5 at different fractions of the viscous timescale t v ). 

From these results we can estimate the true level of angular 
momentum transport present in our disc simulations. We see from 
Table[TJthat for both the 10 5 - and 10 6 -particle runs, there is essen- 
tially no difference in the measured viscosity between the runs with 
Vj n = (i.e., artificial viscosity only) and v- m = 1CT 5 , and that the 
viscosity in the higher-resolution runs is smaller than that in the 
lower-resolution runs (by a factor of approximately 10 1/3 - 2.15), 
suggesting that in this regime artificial transport of angular momen- 
tum is dominant. For larger values of v in , however, the angular mo- 
mentum transport increases as expected, showing that our imposed 
Navier-Stokes viscosity is the dominant source of angular momen- 
tum transport for v in > 10~ 5 (or, equivalently, a > 0.008). An input 
a of 0.01 gives a value of v at R of 1.5 x 10 -4 . The SPH smoothing 
lengths throughout our simulated discs are comparable to those in 
our 10 6 -particle spreading rings at radius R , being of order 0.01 
in code units, and so we use this set of rings for comparison. This 
suggests that the artificial viscosity sets a floor to the effective vis- 
cosity in the SPH simulations, approximately at or slightly below 
our canonical imposed value of a = 0.01. We are therefore satis- 
fied that artificial transport of angular momentum does not domi- 



Table 1. Summary of ring spreading tests. V[ n denotes the magnitude of the 
imposed Navier-Stokes viscosity, as specified in equations l 10l and [T2l while 
v ol ,( is the measured viscosity in the SPH runs, calculated by fitting the 
viscous time-scale as in Equation [T3] The effective a values are calculated 
by assuming that v oul = aCsoHg. For very small values of Vj n the artificial 
viscosity is the dominant source of angular momentum transport, but for 
v in <: 10 5 the measured viscosity increases as expected. 



N Vj n v oll t Corresponding a 



10 5 





3.01 X 10~ 4 


0.019 


10 s 


10- 5 


3.06 x 10~ 4 


0.020 


HP 


10"* 


3.45 x 10~ 4 


0.022 


L0 5 


10- 3 


6.92 x 10~ 4 


0.045 


LO 6 





1.29 x 10~ 4 


0.008 


to 6 


10- 5 


1.38 x 10~ 4 


0.009 


LO 6 


10- 4 


2.25 x 10- 4 


0.015 


LO 6 


10- 3 


1.05 x 10~ 3 


0.068 



nate the viscosity in our disc modelj^. Moreover, it is known that 
in the case of a shearing disc, the SPH arti ficial viscosity behaves 
similarly to a Shakura-Sunyaev a viscosity jMurrav|[l996t) . Conse- 
quently, although the SPH artificial viscosity prevents us from run- 
ning simulations with very low disc viscosities, we are confident 
that numerical angular momentum transport does not dominate our 
results. 



2.2 Numerical set-up and initial conditions 

We choose a system of units (mass Mo, distance ^0 and time Pq) 
such that for a planet mass M P and a stellar mass M„ M P + = 
1M . The unit of time P is the Keplerian orbital period for a semi- 
major axis ^o- This choice of units fixes the gravitational constant 



3 The exception is run PNM, which uses a much lower explicit viscosity. 
In this case we expect the artificial viscosity to dominate the angular mo- 
mentum transport. 
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Table 2. Summary of parameters used in our disc models. 



E [Code Units]^ 


So [g/cm 2 ]^ 


7 


1 


Model Name 


1.12 X 1(T 5 


10 2 


1 


0.005 


Low5 


1.10 X ur 5 


10 2 


1 


0.025 


Low25 


1.12 X ur 4 


io- 5 


1 


0.005 


High5 


i.io x ur 4 


to 3 


1 


0.025 


High25 


i.io x ur 4 


10 3 





0.025 


Flat 


7.03 X 10~ 4 


6.4 x 10 3 





0.025 


PNM^_ 


7.03 x ur 4 


6.4 x 10 3 


1 


0.025 


PNMslope 



a As the unit of mass depends on the planet mass Mp, different values of q 
give a different Xg for the same physical model. 

^ These values correspond to a 1M Q star and an initial semi-major axis of 
1 au. 



c Corresponds to the disc model used in lPNMOll . For this model the vis- 
cosity v in equation [12] was 1.59 x 10~ 6 in dimensionless units, far less than 
the artificial viscosity, see sections ^. 1 . 1 1 & |3. 1 .21 This was set i n error, a nd 
should have been 2.68 X 10~ 5 in our units to match that used in iPNMOll 

to be G = An 2 . This system is particularly convenient as for a mass 
M = IM Q and radius R = lau it gives P = 1 year. 

Our initial conditions consist of a gas disc which is axisym- 
metric about the centre of mass and extends radially from 0.4 to 6 
R . It has a power-law surface density such that 



m) = So 



— 

[Ro 



(14) 



where Eo = E(i?o) is a reference surface density used for normali- 
sation. The viscosity v in equation [T2] for each disc model (barring 
that used in run PNM, see tableO was chosen such that vE was con- 
stant (i.e., the disc is a steady-state accretion disc), and normalised 
such that v = 0.0lc s0 H (wher e the subscript indicates v alues at 
Ro). This treatment reduces to a lShakura & Sunvaevl alpha- 
prescription, with a = 0.01, in the canonical case of a I oc R 1 
surface density profile. The vertical scale-height H is determined 
by the T oc Rr 1 ^ 2 temperature profile described above, so our disc 
has H/R oc R 5/4 . We normalise the temperature profile by setting 
H/R = 0.05 at R = R . 

In each case our disc is modelled using N = 10 7 SPH parti- 
cles. We build the initial conditions by randomly distributing the 
particles in the radial and azimuthal directions according to equa- 
tion[l4] Particles were given zero velocity in the radial and vertical 
(z) directions, and azimuthal velocities such that 



R 



1 dP 
R 2 + p~dk~ 



(15) 



(i.e., Keplerian orbital velocities, with a small correction to account 
for radial pressure gradients). The particles were distributed in the 
vertical direction by randomly sampling a Gaussian density pro- 
file with scale-height H . This set-up is not strictly in equilibrium, 
due to the discontinuities at the inner and outer disc edges, and 
consequently transient density waves pass through the disc for ap- 
proximately one outer orbital period (~ 15J°o)- These transients are 
short-lived, however, and tests performed with an initially relaxed 
disc show that these do not affect our results in any significant way. 

We ran 7 models, using the disc parameters described in table 
|2]and using planet masses q = Mp/M* = 0.005 & 0.025. The star 
and planet were initially placed on the x-axis and given Keplerian 
orbits around the centre of mass, and each model was allowed to 
evolve for 200 planetary orbital periods (unless otherwise stated). 



Artificial Viscosity Only 

Artificial + Navier- Stokes Viscosity 




50 100 

Time (P Q ) 

Figure 2. Comparison between initially eccentric models with just artificial 
viscosity (black line) and with both artificial viscosity and a Navier-Stokes 
viscosity (red line). Both used disc Low25 (see table[2} and were given an 
initial eccentricity eg = 0.05. The model without the Navier-Stokes vis- 
cosity shows small initial damping of eccentricity which soon flattens off, 
while the model with the full viscosity scheme implemented sees continued 
and increased eccentricity damping as it evolves. This shows that the SPH 
artificial viscosity is not causing spurious eccentricity damping. In the case 
of the full Navier-Stokes viscosity model, at later times the eccentricity de- 
cay reversed and began to grow again. This is due to the eccentric planet 
causing stronge r disc e ccentricity, and is in agreement with the findings of 
iD'AngeloetalJbOOfj) . 



The simulations were run on the alice HPC cluster at the University 
of Leicestefl using 128 parallel cores. 



3 RESULTS 
3.1 Code tests 

3.1.1 Numerical eccentricity damping 

As an initial test, and to ensure that the effects we see are physical, 
we first verified that a planet on an eccentric orbit does not undergo 
spurious eccentricity damping due to the SPH artificial viscosity (or 
other numerical effects). To this end we ran 2 realisations the disc 
model Low25, where in each case the planet was given an initial 
eccentricity e = 0.05. In one version of this model the Navier- 
Stokes prescription described in equations[8]to[T0]was switched off, 
so that in this case the only source of angular momentum transport 
was from the artificial viscosity. These were allowed to evolve for 
125 orbits, and the eccentricity evolution is shown in figure[2] 

In the model with no physical (Navier-Stokes) viscosity we 
see very little initial eccentricity decay during the initial period 
while the disc settles into an equilibrium state and the planet opens 
its full gap. A loss of some eccentricity during thi s phase is fully 
expect ed and is in agreement with simulations by iBitsch & Klevl 
who found that for non-gap opening planets eccentricity 
is damped by the surrounding gas. With the physical viscosity 
switched on we see additional damping of eccentricity, at a much 
more pronounced level. In both cases we see exponential damp- 
ing after the initial phase (beyond ~ 75 orbits), after the planet has 
fully cleared its gap. The rate of eccentricity damping during this 



See http : //go . le . ac . uk/alice 
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Figure 3. Surface density evolution of the central region of our disc model PNM (see table[2), roughly equivalent to run N4 from lPNMOll The distance unit 
is equal to the initial separation between the star and planet. Times shown are in units of the initial orbital period of the planet. The eccentricity evolution for 
this model is shown in figure|4] After the inner disc clears (upper panels), the presence of the planet drives the inner edge of the disc eccentric (bottom panels). 
As the disc evolves it exerts torques back upon the planet, causing the planet's eccentricity to grow. [All of the SPH maps presented here were rendered using 
splash lPricell2007t) .l 



phase is similar between the two models. This is because once the 
gap has fully formed, the planet is not directly interacting with the 
gas to any great extent so the angular momentum exchange here 
is due to gravitational resonances. At later times, the eccentricity 
began to ris e again in the case of th e full viscosity model. This is 
expected, as iD'Angelo et alj J2006h found that even for very low 
planet masses an initially eccentric planet can undergo far stronger 
eccentricity growth than one on an initially circular orbit. Conse- 
quently we conclude that the SPH artificial viscosity is not causing 
significant spurious eccentricity damping in our disc models. 



as possible in form to that used in their calculations^. Using the pa- 
rameters for run PNM given in table [2] this approximates run N4 
from that paper, with the obvious caveat that our simulations are in 
3-D. Unlike our other models, the normalisation for the viscosity 
law described in section [2T21 was taken to be 1.59 x 10~ 6 in our di- 
mensionless code units0. This is constant across the disc, fulfilling 
the steady-state accretion requirement that vL be constant. Note, 
however, that with this setup the angular momentum transport due 
to the explicit viscosity is smaller than that due to the SPH artificial 
viscosity (see Section [2.1. IK so in practice our test calculation is 
somewhat more viscous than that of PNM01 (by a factor of a 2-3). 
This model was allowed to evolve for 340 orbital periods. The evo- 



3.1.2 \PNM0 A Result 

A s a furthe r test, we have also attempted to reproduce the results 
of lPNMOll . To this end we have created a disc model that is as near 



5 IPNMOll used a 2-D fixed-grid code for their simulations, so it is not pos- 
sible for us to run a completely identical si mulation . 

6 Equivalent to 2.50 X 1(T 6 in the units of lPNMOll 
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Figure 4. Evolution of planet eccentricity f or the dis c model PNM, which 
is approximately equivalent to ran N4 from PNMOl] (see table [2). We find 
growth of eccentricity in general agreement with that paper. Surface density 
plots from this run are shown in figure [5] The ~ 100 orbital period oscilla- 
tions are due to the relative precession of the planet and the eccentric inner 
edge of the disc. 



lution of the planet's eccentricity is shown in figured while a series 
of surface density maps as the disc evolves are shown in figure [3] 

The evolution of our disc structure is broadly in line with that 
found by IpnmoiL with the planet rapidly opening a wide gap in 
the disc, and the inner part of the disc quickly accreting onto the 
central star. As the system evolves the planet begins to drive ec- 
centricity in the disc at its inner edge, while its own orbit remains 
essentially circular. At later times this is no longer the case and the 
planet's orbit becomes significantly eccentric [ above the ~ 0.01- 
0.05 level required by lOgilvie & Lubowl J2003h for non-coorbital 
corotation resonances to saturate, at which point further eccentric- 
ity growth is expected]. The long-period oscillations in eccentricity 
seen in figure|4]are due to the relative precession of the planet and 
the eccentric disc inner edge of the disc. 

The level of eccentricity gr owth see n in our simulations is 
somewhat less than that seen by IpnmoiL but it is broadly com- 
parable. Moreover, given our larger effective disc viscosity, and the 
inherent differences between the methods (2-D fixed-grid versus 
our full 3-D SPH calculation), we do not expect exact agreement. 
We also note that our simulations are extremely computationally 
expensive (using up to approximately 150,000 CPU hours per run), 
so we are limited in how long we can evolve our models for. Conse- 
quently we are unable to find a level of eccentricity at which growth 
saturates (the eccentricity was still growing at the end of our sim- 
ulation), but otherwise we find g ood agre ement between our 3-D 
results and the 2-D simulations o f lPNMOll . 



3.2 The effect of surface density & planet mass 

3.2.1 Planet mass 

To test the effect of different planet masses, we ran models with 
both Low and High surface densities (see table [2} with planet-star 
mass ratios of q = 0.005 & 0.025. For a 1M G star this corresponds 
to planet masses of approximately 5 and 25Mj up respectively. Fig- 
ure [5] shows the time evolution of the orbital eccentricity in each 
case. For both of the models with q = 0.005 (Low5 and High5), 
no eccentricity growth was seen at any level. For the models with 




200 



Figure 5. Eccentricity evolution for planets of different masses in disc mod- 
els with different surface densities. Low5 and High5 have a planet-star mass 
ratio of q = 0.005, while Low25 and High25 have q = 0.025. Low and High 
refer to the choice of disc surface density - see table [2]for the values - but 
all have a power-law index of y = 1 . 



q = 0.025 (Low25 and High25) some modest growth was seen, but 
none above e = 0.005. 

This is consistent with the conclusions of both IpnmoiI and 
iD'Angelo et aTlfeOOfjl) . who found that eccentricity is first induced 
in the disc, driven by the outer 3:1 Lindblad resonance. We clearly 
see this disc eccentricity in both runs with q = 0.025, but not in 
either of the lower mass planet cases. A comparison of the surface 
density maps for both planet masses in the Low model is shown 
in figure [6] The reason for the lack of eccentricity growth in both 
of the q = 0.005 models is that the 3:1 outer Lindblad resonance 
induced by the planet is too weak to affect the disc structure sig- 
nificantly. The difference in eccentricity evolution between models 
Low25 and High25 arises because the higher surface density disc 
is simply more massive, and consequently is able to exert stronger 
torques upon the planet, resulting in more (but still extremely lim- 
ited) eccentricity growth. 

We are confident that in the case of both planet masses, eccen- 
tric co-rotation resonances are resolved where prese nt. Using the 
width formula provided bv lMasset & Ogilviej ( 120041) . our smooth- 
ing lengths are smaller than the resonant widths at the nominal res- 
onant locations by a factor of at least two, for resonances that are 
not fully in the open gap. For the case of the 5 M Jup planet, this is 
assuming an eccentricity of 10~ 2 , as for the eccentricity measured 
from the simulations the prescribed widths are vanishing. 



3.2.2 Radial surface density profile 

We now consider the effect of the disc surface density profile on the 
evolution of the embedded planet. To this end we have run two ad- 
ditional models with the same star-planet mass ratio of q = 0.025, 
and different radial surface density profiles (y = 1 & 0). The 
disc models were Flat and PNMslope (see tabled, and each was 
evolved for 200 orbits. The evolution of the orbital eccentricity for 
these runs, along with two previously described (High25 and PNM, 
for the purposes of comparison), are shown in figure|7j Again, due 
to computational limitations we have not run these models for long 
enough to determine at what level of eccentricity growth saturates; 
we are instead more concerned here with determining the condi- 
tions under which growth will occur. We see again the eccentric in- 
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Figure 6. Surface density maps after 200 planetary orbits: panel (a) shows 
the central region of model Low5 (q = 0.005) and panel (b) shows Low25 
(q = 0.025; see table|2j. It is clear from panel (a) that the lower mass planet 
has a far smaller perturbing effect on its host disc, driving spiral waves but 
not further disrupting the disc shape or structure. This is in stark contrast to 
panel (b), where the higher mass planet has driven the inner edge of the disc 
quite eccentric. 



ner disc driven by the presence of the massive planet, as described 
above, and again find that its effect on the eccentricity of the planet 
depends strongly on the precise disc model used. 

The results from these models show two things. First, the mag- 
nitude of the disc surface density has a strong effect on the eccen- 
tricity evolution of an embedded planet. Figure [7] clearly shows 
that the two models with S = 7.03 x 1CT 4 (PNM and PNMs- 
lope) see significant eccentricity growth, while the models with 
So = 1-10 x 10~ 4 (High25 and Flat) do not. This suggests that 
eccentricity can only be excited a bove a thr eshold surface density. 
This behaviour was suggested bv lPNMOll but not investigated in 
detail. Comparing the ratio between the planet mass and the local 
disc mass (approximated by ~Ljxa 2 in the unperturbed disc) suggests 
that values between ~ 1 - 13.5 may result in eccentricity growth 
(see figure [SJ. We also note in passing that this is a surprisingly 
strong effect for a relatively modest (factor of 6.4) change in the 
disc surface density; the surface densities in real protoplanetary 
discs are thought to change by factors > 10 3 over their lifetim es 
(e.g jHartmann et al.lll998l ; rAlexander. Clarke. & Pringlej|2006h . 

The power-law index y also has quite a strong effect on the 
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Figure 7. Eccentricity evolution for disc models with different surface den- 
sity profiles. The black line is the same as in figure|4] truncated at t = 200Po 
for reference, and the green line is the same as in figure \5\ These models 
show that while higher values of give consistently stronger eccentricity 
growth, shallower radial profiles (lower values of y, see table |5J also show 
stronger growth. The oscillations in the curves are again due to the rela- 
tive precession of the planet and the eccentric disc. Note also that we ran 
model PNMslope (red line) for a further 100 orbits (not shown here) to ver- 
ify that the eccentricity does grow indeed continue to grow; the declining 
eccentricity at t = 200 is simply due to this precession effect. 



level of eccentricity growth seen. The two models with y = 1 
(High25 and PNMslope) show slower, weaker growth of eccentric- 
ity than their counterparts with the same normalisation value of Eo 
but a flat 7 = radial power-law dependance. There are two pri- 
mary reasons for this. Firstly, a flatter radial profile puts less mass 
inside the planet's orbit. The inner disc therefore accretes on to the 
star more rapidly for flatter surface density profiles (i.e., lower val- 
ues of y), and as eccentricity only begins to grow considerably after 
the inner disc has accreted, this takes place sooner for flatter surface 
density profiles. A flatter radial profile also puts more mass into the 
outer Lindblad resonances, which are responsible for the torques 
that drive eccentricity growth. This changes the torque balance on 
the planet, and gives rise to more rapid eccentricity growth. In addi- 
tio n, in the special case wh en Lindblad torques cancel (as proposed 
bv lGoldreich & Saril2003l) . the resultant net torque is a function of 
the surface density gradient, with a steeper power-law dependance 
damping eccentricity. However, we do not expect this effect to be 
active in our models, as the planets open gaps sufficiently wide that 
co-rotation resonances are ineffective. 

To investigate this effe ct further we have followed the method 
of lArtvmowicz et all J 199 lh to calculate values of e, time averaged 
over several orbits at the end of each simulation. The radial con- 
tributions to this, normalised against the magnitude of the surface 
density in each model, are shown in figure [9] for all models with 
q = 0.025. In the low-surface density limit (models Low25, High25 
and Flat), we see that the eccentricity is being damped by a peak 
at R ~ 1.8, and the magnitude of the surface density acts as a scal- 
ing factor with only a very weak dependence on the radial slope 
(y). In these three disc cases, e is always negative. By contrast, 
In the models with higher surface densities, where eccentricity is 
growing (models PNM and PNMslope), the sense of the contribu- 
tion reverses, and the total e is positive. Here there is also a more 
pronounced effect of the radial gradient in the surface density. It is 
unclear exactly what mechanism causes this effect, but it appears 
that there is a threshold surface density above which the analysis 
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Figure 8. Comparison between the local unperturbed disc mass (approxi- 
mated by ixZc?) and planet mass for simulations presented here and in two 
other pa pers. Crosses indic ate simulations f r om th is paper, squares from 
IPNMOll and triangles from 15] Angelo et "all <2006h . Green symbols indi- 
cate that eccentricity growth was seen, while red indicates that it was not. 
The choice of units on the x-axis gives the reference surface density at 
the semimajor-axis of the planet for that model (except for the models of 
iD'Angelo et alj fefJOd) . which placed their planets at 5.2au; in this case the 
x-axis value is 5.2 2 times their reference surface density). The dashed lines 
show where the ratio of M p /nl.a 2 cross values of 1.0 and 13.5. We sug- 
gest that these values represent cases where the planet is massive enough to 
significantly perturb the inner edge of the disc, but where the disc is also 
massive enough that it can exert sufficiently strong torques upon the planet, 
although we note that this ratio is not the lone deciding factor. We also draw 
attention to the fact that the single point below the threshold M„lKLct = 1.0 
has been questioned bv lMasset & Qgilviel j2004h as being affected by spu- 
rious numerical factors, and thus the lower limit may not be real. We do not 
expect this to extend into the sub-jovian regime, where gap-opening is not 
necessarily efficient enough to allow similar modes of eccentricity growth 
(although this also depends on other factors including disc viscosity). 



perfo rmed in the low disc mass limit (e.g. iGoldreich & Tremainel 
Il98d) no longer applies. Our results, and those of previous simula- 
tions (see Figure[8}, suggest that this threshold is given by 

nla 2 > M p /C (16) 

where C is a constant with a value C ~ 10. Discs with surface den- 
sities below this threshold are unable to excite significant eccentric- 
ity in the planet's orbit. This behaviour was suggested by PNM01, 
and our results support their tentative prediction. The threshold 
must presumably depend on several other factors (disc viscosity, 
H/R, etc.), and a complete exploration of this parameter space is 
beyond the scope of this work. Nevertheless, our results show that 
eccentricity is only excited in discs with very high surface densities. 



4 DISCUSSION 

4.1 Numerical limitations 

The biggest potential numerical problem with this work is that the 
SPH artificial viscosity may cause spurious eccentricity damping. 
Artificial viscosity can be especially problematic for shearing-disc 
type problems such as this, as the differential rotation is often mis- 
taken by the algorithm for a shock. The effects of artificial viscosity 
typically scale with the SPH smoothing length h, and consequently 
are a strong function of numerical resolution. However, we have 



been able to run our simulations at very high resolution, and given 
the results of the various tests presented in sections l2. 1 . 1 I andl 3 . 1.11 
we are confident that the SPH artificial viscosity is not a dominant 
influence on our results. 

We further note that our approximation of a locally isothermal 
equation of state is somewhat idealised, and in particular may not 
give an accurate descrip t ion of the spiral density waves induced in 
the disc. iBitsch & Klevl fcOlCh looked at the evolution of initially 
eccentric planets in fully radiative discs, and made comparisons to 
models which used an isothermal approximation. They only found 
the results to be inconsistent for relatively low planet masses (< 
0.6 Mj up ). The planet masses considered here are far above this, 
well into the regime where the isothermal approximation holds, and 
so we do not expect our use of an isothermal equation of state to 
introduce significant uncertainties in our results. 

Perhaps more of a concern is the viscosity in our simulated 
discs. We are confident that the SPH artificial viscosity in our sim- 
ulations does not dominate our results, but it does set a floor to 
the range of viscosities we can explore: with current computa- 
tional capabilities we cannot model discs with a < 0.0 fl Pre- 
vious simulations using grid-based methods have adopted a wide 
range of viscosities and temperature prescr iptions, with effective a 
values ranging from a ~ 10 4 — 10 2 (e.g., Papaloizou et al.ll200ll : 



ID'Angelo et"al1l2006t IBitsch & Kley||2010l) . It is not clear which 
value(s) of a is most appropriate, but we note that our simula- 
tions have a viscosity that is somewhat larger (by a factor = 2-3) 
than the canonical values adopted in previous (mostly 2-D) studies. 
We note also that values inferred from observatio ns of protoplan- 
etary discs are typically of order a ~ 0.01 (e.g. lHartmann et al .1 
1 19981 : iKing. Pringle. & Lividl2007t) . Presumably the surface den- 
sity threshold for eccentricity growth must depend on the disc vis- 
cosity (and indeed temperature), but unfortunately we are not able 
to explore this issue further. 

It must also be borne in mind that the Navier-Stokes viscos- 
ity used in our models is merely a first-order approximation, at- 
tempting to mimic the effect globally of a process that occurs far 
below the scales we are able to resolve - namely the magnetorota- 
tional instability (MRI) thoug ht to drive angular mom entum trans- 
port in protoplanetary discs jBalbus & Hawlev| [l991). By adopt- 
ing a lShakura & Sunvaevi d 19731) alpha-prescription we implicitly 
assume that the small-scale effects of turbulence in the disc be- 
have like a viscosity on large scales, but it is not clear whether 
this approximation holds for length scales < H. In our simula- 
tions the planets open gaps on length scales > H, so turbulent 
fluctuations on smaller scales are unlikely to have a strong effect. 
However, there is some overlap between the length-scales consid- 
ered here and the typical scales of MRI turbulence, so we note 
that our results may not hold if the MRI drives significant power 
on moderate or large scales (as sug gested by recent simulations; 
ISimon. Beckwith. & Armitagel2012l) . Detailed investigation of this 
issue is beyond the scope of this paper, but if the viscous approxi- 
mation does break down at scales ~ H then it seems likely that ec- 
centricity growth could be affected, particularly for low-mass plan- 
ets. 



7 Note that angular momentum tr ansport by the SPH artificial viscosity 
scales with the smoothing length /; (Murray 199oJ), which in 3-D scales as 
h oc A? 1 ' 3 . Reducing the numerical viscosity to a < 0.001 would therefore 
require us to increase the SPH particle number N by a factor ~ 10 3 , which 
is not currently feasible. 
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Figure 9. Radial contributions to e, c alculated using a Gaussia n perturba- 
tion approach following the method of Artymowicz et al. ( 1991). These val- 
ues are time-averaged over 5 orbital periods of the planet at the end of each 
simulation. Models Low25, High25 and Flat clearly show that in the low- 
surface density limit, e is approximately linear with the magnitude of the 
surface density, with a weak dependance on the radial profile y. For models 
PNM and PNMslope it can be seen that above the threshold for eccentricity 
growth, the sign of the curves flip and the net e becomes positive. In this 
limit the radial profile of the surface density becomes more important, and 
the effect is no longer linear with its magnitude. 



4.2 Applications to real systems 

Our major result is that resonant torques are not generally an ef- 
ficient means of exciting eccentricity, and that planet-disc interac- 
tions are unlikely to be responsible for the eccentricities seen in 
the majority of exoplanet systems. While we do find eccentricity 
growth in agreement with lPNMOll we note that their calculations 
considered a very massive planet in a disc with a very high surface 
density. Their reference surface density of 6.4 x 10 4 g/cm 2 at 1 au is 
a factor of several larger than pred icted by the Minimum Mass So- 
lar Nebula dWeidenschillin dl 1 9771) or more realistic accretion disc 

We find that a modest reduc- 



models (e.g.. lHartmann et al 



tion in the disc surface density results in no significant eccentricity 
growth for si milarly m assive planets, in agreement with behaviour 
suggested bv lPNMOll Consequently it is unlikely that this mecha- 
nism will be able to excite eccentricity in real protoplanetary discs. 

We also fail to find eccentricity growth for lower planet 
masses, but in this cas e we treat our results with more caution. 
Using 2-D simulations iD'Angelo et al.l d2006h found eccentricity 
growth up to e ~ 0.1 for planets between 1-3 Mj up , but typically 
this occurred on time-scales of thousands of orbits. Given the high 
computational cost of our 3-D simulations we are not able to follow 
their evolution for such long time-scales, and consequently we can- 
not rule out this slo wer mode of growth. Ho wever, we note that the 
disc model used by ID'Angelo et alj ( 120061 ) also uses a very large 
surface density: if we re-scale their model to match the units used 
here, their surface density normalisation (E ) becomes 2.0 x 10 3 
g/cm 2 at 1 au, larger than in our High models. We also note that 
their relatively flat choice surface density profile (y = 1 /2) is likely 
to promote eccentricity growth. Although our simulations do not 
allow us to rule out growth on very long time-scales, analysis like 
that shown in figure [9] shows that the disc gives a negative con- 
tribution to the planet's e which indicates that eccentricity will 
not grow unless the disc structure changes significantly. Instead it 
seems likely that differences in the choice of disc model are respon- 



sible for the apparent d iscrepancy between our results and those of 
ID'Angelo etail feOOfjl) . 

As noted in section [3~.2.2| there are two reasons why the radial 
surface density profile plays a role in exciting eccentricity growth. 
At a basic level, a lower value of y puts mass mass into the outer 
Lindblad resonances and conversely less mass into the inner reso- 
nances. A s the mode of growth relies on the 1:3 outer resonance 
dPNMOll) . this favours eccentricity growth by strengthening that 
resonance. In the case of lower mass planets, where a gap is not 
fully opened, there is another effect which takes place that brings 
the ra dial surface density gradient into play. iGoldreich & S aril 
d2003l) suggest that in a near-Keplerian disc where an outer and in- 
ner Lindblad resonance cancel to a reasonable approximation, the 
resulting net torque scales with dZ/dr rather than with E. This im- 
plies that not only the S level but also the radial profile may be- 
come as important as the planet mass in discerning the physical 
contribution of resonant torques for low mass planets. We believe 
that it is the former effect that we are seeing in our simulations, 
with both the surface density profile and its normalisation level both 
having a strong effect on if, and how, the orbital eccentricity of an 
embedded planet will grow (figure[7](. 

This analysis of the effect of both the magnitude and radial 
profile of the surface density is supported by looking at the ra- 
dial contributions to the change in eccentricity in different models, 
shown in figure [9] Below the threshold for eccentricity growth the 
surface density profile has only a small effect, and its magnitude 
is an approximately linear scaling factor, wh ich is in agreement 
with the findings of I Artymowicz et al] J 1 99 lh . In contrast, above 
the threshold for growth the sense of the contributions reverse and 
the difference between models with the same So but different y be- 
comes more pronounced. We are confident that this effect is real, 
and warrants further study, but detailed investigation is beyond the 
scope of this paper. Other disc parameters including viscosity and 
temperature structure must also have an effect on this result, but we 
have not investigated them here. 

We can extend this analysis by taking equation[T6]to be a nec- 
essary condition for eccentricity growth: 



CnLa- > Mn 



(17) 



where C ~ 10 is a constant. However, giant planets continue to 
accrete via tidal streams even after opening a gap in the disc, and 
in discs which are sufficiently massive to meet this condition they 
are likely to accrete rapidly. This accretion increases M p , making 
it less likely that the threshold for eccentricity growth will be met. 
The critical issue is therefore one of time-scales: does eccentricity 
grow on a shorter time-scale than the planet's mass? For massive 
giant planets (> 5Mj up ) tidal torq ues strongly suppress accretion 
from the disc on to the planet (e.g. jLubow. Seibert. & Artvmowiczl 
1 19991) , and we therefore expect eccentricity to grow more rapidly 
than the planet mass. However, our results suggest that eccentricity 
growth in this regime requires very high disc surface densities, > 
10 3 gctrT 2 (see Fig[8}, and such massive discs are not commonly 
observed (see below). 

By contrast, lower-mass planets (~ 1M Jud ) accrete very effi- 
cientl y from their parent discs (e.g.. ID'Angelo. Henning. & Klevl 
|2002|) . We can parametrize the planetary accretion rate as M p = 
fM d , where M A = 3ttvE is the disc accretion rate in the absence of 
the planet. The parameter e represents the efficiency of accretion on 
to the planet; simulations show that this efficiency has a peak value 
of e ~ 1 for planets o f approximately Jupiter mass, and decline s to 
higher planet masses jLubow et alj 1999| - |D'Angelo et al.l2002l see 
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also lVeras & Armitagell2004l) . If we substitute v = aQ.H- 
rearrange equation [T6] to find 



Mp < C -io- 
—r- < — a 12 

M n 3e 



(!)" 



(18) 



The quantity on the left-hand-side is the time-scale for planet 
growth through accretion of gas from the disc, T acclctc , so we see 
that meeting the condition for eccentricity growth (equation I16t 
also sets an upper limit to the time-scale for planet growth by ac- 
cretion. Taking standard values of a = 0.01 and H/a = 0.1, and 
assuming that e ~ 1, we find 



f accrete 5; 3 X 10 £2 , 



(19) 



for planets of approximately Jupiter mass. This is comparable to 
the time-scale for ec centricity growth seen in previous studies of 
Jupiter-mass planets dD' Angelo et alj|2006l) . We therefore suggest 
that in this case the planet would accrete rapidly, and 'migrate' out 
of the region of allowed eccentricity growth highlighted in figure [8] 
before it attains a significant eccentricity. 

We have shown that in moderately viscous discs (with a ~ 
0.01), eccentricity growth only occurs if the disc surface density is 
high. Unfortunately, observational constraints on the surface densi- 
ties of discs are extremely weak, and at au scales such as those con- 
sidered here, almost non-ex istent. [Andrews & Williams] d2007l) and 
lAndrews et alj |2009, 2010) made systematic studies of protoplan- 
etary discs in Taurus and Ophiuchus, fitting surface density profiles 
to submillimeter observations of thermal dust emission. However, 
the angular resolution of such observations means that they can 
only resolve scales a few tens of au (in the very best cases), and are 
not sensitive to the region of most interest for planet formation. If 
we extrapolate their results down to the unresolved inner disc, their 
fits yield values of So at au radii between ~ 30 and ~ 1700 g/cm 2 , 
with most values lying at ~ 700 - 800 g/cm 2 . The radial power- 
law indices (y) they fit to their observations range between 0.4 and 
1.1, with a clear preference towards the upper end of this range. 
These correspond approximately to our High disc model, and have 
both lower surface d ensities a nd s teeper power-law indic es than the 
disc models used bv lPNMoil and lD' Angelo et alj bOOd) . However 
it must be stressed that extrapolating these sub-mm observations 
is by no means justified. It has even been suggested that substan- 
tial mass reservoirs may exist in 'dead zones' close to the star (e.g., 
lGammiell99l ; lHartmann et aLfeOOelfzhu et al.l2010t) . but there are 
currently no useful constraints on disc surface densities at au radii. 
ALMA may soon provide real constraints on protoplanetary discs 
with far higher angular resolution, but it will be some time before 
it is a ble to probe radii of a few au (e.g.. lCossins. Lodato. & TestH 

[2oToh . 

Despite this, we argue that significant eccentricity growth due 
to the planet-disc interaction is unlikely given realistic protoplane- 
tary disc conditions. There is an increasing consensus i n the litera- 
ture that this is the case - while lD' Angelo et alj | |2006|) did see ec- 
centricity growth by this method, it did not rise above ~ 0. 15, lower 
than the observed values of 0.2-0.3 that thi s effect is invoked to 
explai n. Similarly, the semi-analytic models rMoorhead & Adami 
J2008I) were unable to reproduce the observed low-eccentricity dis- 
tribution. Coupled with our findings, it seems that the planet-disc 
interaction alone is incapable of reproducing the eccentricities seen 
in exoplanet observations. 

This negative result of course begs the question of what the 
true origin of exoplanet eccentricities is. An emerging consensus 
seems to be that scatteri ng events are responsib le for most if not all 
of the distribution (e.g. IChatteriee et alj|2008l) . This is backed up 



by a vast n umber of tightly-pa cked multi-planet systems observed 
by Kepler teatalha et al]|20i2l) and by a plethora of planetary sys- 
tems that have clearly undergo ne strong interactio ns (e.g. highly in- 
clined and retrograde planets; IWright et al .11201 j) . However, there 
is some suggestion that the low-eccentricity end of the distribution 
perhaps cannot be explained b y this mechanism jGoldreich & Saril 
l2003l : |juric & Tremainell2008h . and further work in this area is still 
required. 



5 CONCLUSIONS 

We have performed high-resolution 3D SPH simulations of giant 
planets embedded in protoplanetary discs. For high disc surface 
densities and planet masses we find that the planet-disc interac- 
tion leads to eccentricity growth, in agreement with previous stud- 
ies. However, we have shown that for realistic planet masses and 
disc properties, the planet-disc interaction is incapable of exciting 
significant orbital eccentricity growth. While we do not investigate 
the effect of different disc viscosities, we identify a threshold sur- 
face density for eccentricity growth, and note that this threshold 
is rarely, if ever, met in real systems, except in cases where the 
timescale for eccentricity growth is comparable to the timescale 
for mass accretion by the planet. We conclude from our simula- 
tions that in the case of a real giant planet, the interaction with its 
parent disc is unlikely to yield growth of its orbital eccentricity at 
measurable levels. Therefore we suggest that the low but non-zero 
exoplanet eccentricities observed, not accounted for by simulations 
of planet-planet scattering events, must have some other origin. 
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